Land Cover 2010 to 2020#

From 2010 to 2020, the Denver MSA has maintained control of urban sprawl. While the population increased by 16%, from 2.54 million to 2.95 million (block-level population from the 2010 and 2020 decennial census datasets), the urban area has expanded by only 3.5%, from 1,200 to 1,240 square miles, according to data from MRLC.

Hide code cell content
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio as rio
from shapely.geometry import Polygon, Point
crs = "EPSG:2232"
%run assets/colors.py
%run assets/rcParams.py
C:\Users\leeje\OneDrive\G-4\land-use-env-modeling\assignments\sprawl-allocation\memo\assets\rcParams.py:8: MatplotlibDeprecationWarning: The seaborn styles shipped by Matplotlib are deprecated since 3.6, as they no longer correspond to the styles shipped by seaborn. However, they will remain available as 'seaborn-v0_8-<style>'. Alternatively, directly use the seaborn API instead.
  plt.style.use("seaborn-whitegrid")

We used land cover data from 2008 and 2019 as proxies of the land cover in 2010 and 2020, as these are the best available data from MRLC. Land cover data rasters with 4000*4000 sqft cell sizes. Our study identified 5 types of land cover types:

  • Urban.

  • Forests.

  • Farms.

  • Wetlands.

  • Water.

Hide code cell content
# Read land cover rasters
# The file paths are relative to the notebook source file's location in the repo's main branch
land_cover_2009 = rio.open("../data/land-cover/2009.tif")
land_cover_2019 = rio.open("../data/land-cover/2019.tif")

def raster_to_shape(raster, crs, to_point=False, value_column="value"):
    data_array = raster.read(1)
    transform = raster.transform
    shapes_data = []
    height, width = data_array.shape

    if to_point:
        for row in range(height):
            for col in range(width):
                value = data_array[row, col]
                if np.isnan(value):
                    continue
                x, y = transform * (col + 0.5, row + 0.5)
                # Create a polygon for the cell
                point = Point(x, y)
                # Add the polygon and its value to the list
                shapes_data.append({"geometry": point, value_column: value})
        return gpd.GeoDataFrame(shapes_data, crs=crs)

    for row in range(height):
        for col in range(width):
            value = data_array[row, col]
            if np.isnan(value):
                continue
            x_min, y_max = transform * (col, row)
            x_max, y_min = transform * (col + 1, row + 1)
            # Create a polygon for the cell
            polygon = Polygon(
                [(x_min, y_min), (x_min, y_max), (x_max, y_max), (x_max, y_min)]
            )
            # Add the polygon and its value to the list
            shapes_data.append({"geometry": polygon, value_column: value})
    return gpd.GeoDataFrame(shapes_data, crs=crs)
Hide code cell content
# Get a fishnet of info from the land cover rasters
fishnet = (
    raster_to_shape(
        land_cover_2009, crs, to_point=False, value_column="land_cover_2009"
    )
    .query("land_cover_2009 != 0")
    .copy()
)
points_2019 = (
    raster_to_shape(land_cover_2019, crs, to_point=True, value_column="land_cover_2019")
    .query("land_cover_2019 != 0")
    .copy()
)
# Join points to fishnet, no index right
fishnet = gpd.sjoin(fishnet, points_2019, how="left", predicate="contains").drop(
    columns=["index_right"]
)

# Classify land cover types
developed = [21, 22, 23, 24]
land_covers = {
    "developed": [21, 22, 23, 24],
    "forest": [41, 42, 43],
    "farm": [81, 82],
    "wetland": [90, 95],
    "water": [11],
}
for year in [2009, 2019]:
    # First update whether developed bool
    fishnet[f"developed_{year}"] = fishnet[f"land_cover_{year}"].isin(developed)

    # Then update land cover
    fishnet[f"land_cover_type_{year}"] = "other"
    for land_cover, codes in land_covers.items():
        fishnet.loc[
            fishnet[f"land_cover_{year}"].isin(codes), f"land_cover_type_{year}"
        ] = land_cover

    fishnet.drop(columns=[f"land_cover_{year}"], inplace=True)
Hide code cell content
import altair as alt
alt.data_transformers.disable_max_rows()

def altair_fishnet(fishnet, column, color_dict, legend_title, title):
    chart = alt.Chart(
        fishnet.to_crs(4326)[["geometry", column]]
    ).mark_geoshape().encode(
        color=alt.Color(
            f"{column}:N",
            title=legend_title,
            scale=alt.Scale(
                domain=list(color_dict.keys()), range=list(color_dict.values())
            ),
        )
    ).properties(
        width=400, height=330, title=title
    )
    return chart
Hide code cell content
from colors import (
    palette_primary,
    palette_backup2_dimmed,
    palette_backup1,
    palette_green,
    palette_green_faded,
    palette_water,
    palette_hero,
)
color_dict = {
    "developed": palette_hero,
    "forest": palette_backup2_dimmed,
    "farm": palette_green,
    "wetland": palette_green_faded,
    "water": palette_water,
    "other": "#dddddd",
}
Hide code cell source
altair_fishnet(
    fishnet=fishnet,
    column="land_cover_type_2009",
    color_dict=color_dict,
    legend_title="Land Cover Type",
    title="Land Cover in 2008",
)
Hide code cell source
altair_fishnet(
    fishnet=fishnet,
    column="land_cover_type_2019",
    color_dict=color_dict,
    legend_title="Land Cover Type",
    title="Land Cover in 2019",
)
Hide code cell source
fishnet["change"] = "Undeveloped"
fishnet.loc[fishnet["developed_2009"], "change"] = "Developed by 2010"
fishnet.loc[
    (fishnet["developed_2019"] & (~fishnet["developed_2009"])), "change"
] = "Developed 2010-2020"

color_dict = {
    "Developed by 2010": palette_hero,
    "Developed 2010-2020": palette_primary,
    "Undeveloped": "#dddddd",
}

altair_fishnet(
    fishnet=fishnet,
    column="change",
    color_dict=color_dict,
    legend_title="Development Change",
    title="Development Change 2010-2020",
)